In [1]:
import glob
import matplotlib.pyplot as plt
import pickle
import lmfit
import numpy as np
import scipy.stats as ss
plt.style.use('seaborn-white')
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from statsmodels.formula.api import ols, rlm
import sys
sys.path.append('../')
from Linearity import Neuron
In [33]:
filename = '/media/sahil/NCBS_Shares_BGStim/patch_data/161013/c1/plots/c1.pkl'
In [34]:
control_result2_rsquared_adj = []
control_result1_rsquared_adj = []
control_var_expected = []
gabazine_result2_rsquared_adj = []
gabazine_result1_rsquared_adj = []
gabazine_var_expected = []
tolerance = 5e-4
In [41]:
def linearModel(x, a=1, b=1):
# Linear model
return (a*x+b)
def DN_model(x,a=1,b=1):
# Divisive normalization model
#return x - a*(x**2)/(b+x)
return x/(a*x+b)
def DN_1_model(x,a=1):
# Divisive normalization model
#return x - a*(x**2)/(b+x)
return a*x/(x+1)
def DN_2_model(x,a=1,b=1):
# Divisive normalization model
#return x - a*(x**2)/(b+x)
return x/(a*x+b)
In [42]:
lin_aic = []
dn_aic = []
lin_chi = []
dn_chi = []
control_observed = {}
control_observed_average = {}
gabazine_observed ={}
gabazine_observed_average = {}
control_expected = {}
control_expected_average = {}
gabazine_expected ={}
gabazine_expected_average = {}
feature = 0
neuron = Neuron.load(filename)
for expt in neuron.experiment:
print ("Starting expt {}".format(expt))
for numSquares in neuron.experiment[expt].keys():
print ("Square {}".format(numSquares))
if not numSquares == 1:
nSquareData = neuron.experiment[expt][numSquares]
if expt == "Control":
coords_C = nSquareData.coordwise
for coord in coords_C:
if feature in coords_C[coord].feature:
control_observed_average.update({coord: coords_C[coord].average_feature[feature]})
control_expected_average.update({coord: coords_C[coord].expected_feature[feature]})
control_observed.update({coord: []})
control_expected.update({coord: []})
for trial in coords_C[coord].trials:
if feature in trial.feature:
control_observed[coord].append(trial.feature[feature])
control_expected[coord].append(coords_C[coord].expected_feature[feature])
elif expt == "GABAzine":
coords_I = nSquareData.coordwise
for coord in coords_I:
if feature in coords_I[coord].feature:
gabazine_observed.update({coord: []})
gabazine_expected.update({coord: []})
gabazine_observed_average.update({coord: coords_I[coord].average_feature[feature]})
gabazine_expected_average.update({coord: coords_I[coord].expected_feature[feature]})
for trial in coords_I[coord].trials:
if feature in trial.feature:
gabazine_observed[coord].append(trial.feature[feature])
gabazine_expected[coord].append(coords_I[coord].expected_feature[feature])
print ("Read {} into variables".format(filename))
In [43]:
list_control_observed = []
list_gabazine_observed = []
list_control_expected = []
list_gabazine_expected = []
if len(gabazine_observed):
for key in gabazine_observed.keys():
for element1, element2 in zip(gabazine_observed[key], gabazine_expected[key] ):
if not (element1 <0 or np.isclose(element1, 0, atol=tolerance) or element2<0 or np.isclose(element2, 0, atol=tolerance)):
list_gabazine_observed.append(element1)
list_gabazine_expected.append(element2)
if len(control_observed):
for key in control_observed.keys():
for element1, element2 in zip(control_observed[key], control_expected[key] ):
if not (element1 <0 or np.isclose(element1, 0, atol=tolerance) or element2<0 or np.isclose(element2, 0, atol=tolerance)):
list_control_observed.append(element1)
list_control_expected.append(element2)
In [44]:
list_control_observed = []
list_gabazine_observed = []
list_control_expected = []
list_gabazine_expected = []
if len(gabazine_observed_average):
for key in gabazine_observed_average.keys():
element1, element2 = gabazine_observed_average[key], gabazine_expected_average[key]
if not (element1 <0 or np.isclose(element1, 0, atol=tolerance) or element2<0 or np.isclose(element2, 0, atol=tolerance)):
list_gabazine_observed.append(element1)
list_gabazine_expected.append(element2)
if len(control_observed_average):
for key in control_observed_average.keys():
element1, element2 = control_observed_average[key], control_expected_average[key]
if not (element1 <0 or np.isclose(element1, 0, atol=tolerance) or element2<0 or np.isclose(element2, 0, atol=tolerance)):
list_control_observed.append(element1)
list_control_expected.append(element2)
In [45]:
minPoints = 10
minIQR = 3
if len(list_control_expected)>minPoints and len(list_control_observed)> minPoints and ss.iqr(list_control_expected)>minIQR:
print ("Control")
X = np.array(list_control_expected)
y = np.array(list_control_observed)
idx = np.argsort(X)
X = X[idx]
y = y[idx]
linear_Model = lmfit.Model(linearModel)
DN_Model = lmfit.Model(DN_model)
lin_pars = linear_Model.make_params()
lin_result = linear_Model.fit(y, lin_pars, x=X)
lin_aic.append(lin_result.aic)
lin_chi.append(lin_result.redchi)
DN_pars = DN_Model.make_params()
DN_result = DN_Model.fit(y, DN_pars, x=X)
dn_aic.append(DN_result.aic)
dn_chi.append(DN_result.redchi)
print (lin_result.fit_report())
print (DN_result.fit_report())
ax = plt.subplot(111)
ax.scatter(X, y, alpha=0.2)
ax.set_xlim(xmin=0.)
ax.set_ylim(ymin=0.)
ax.set_xlabel("Expected")
ax.set_ylabel("Observed")
ax.set_title("Divisive Normalization and Inhibition fits")
ax.plot(X, lin_result.best_fit, '-', label="Divisive Inhibition Model")
ax.plot(X, DN_result.best_fit, '-', label="Divisive Normalization Model")
plt.legend()
plt.show()
In [46]:
if len(list_gabazine_expected)>minPoints and len(list_gabazine_observed)>minPoints and ss.iqr(list_gabazine_expected)>minIQR :
print ("GABAzine")
X = np.array(list_gabazine_expected)
y = np.array(list_gabazine_observed)
idx = np.argsort(X)
X = X[idx]
y = y[idx]
linear_Model = lmfit.Model(linearModel)
DN_Model = lmfit.Model(DN_model)
lin_pars = linear_Model.make_params()
lin_result = linear_Model.fit(y, lin_pars, x=X)
lin_aic.append(lin_result.aic)
lin_chi.append(lin_result.redchi)
DN_pars = DN_Model.make_params()
DN_result = DN_Model.fit(y, DN_pars, x=X)
dn_aic.append(DN_result.aic)
dn_chi.append(DN_result.redchi)
print (lin_result.fit_report())
print (DN_result.fit_report())
ax = plt.subplot(111)
ax.scatter(X, y, alpha=0.2)
ax.set_xlim(xmin=0.)
ax.set_ylim(ymin=0.)
ax.set_xlabel("Expected")
ax.set_ylabel("Observed")
ax.set_title("Divisive Normalization and Inhibition fits")
ax.plot(X, lin_result.best_fit, '-', label="Divisive Inhibition Model")
ax.plot(X, DN_result.best_fit, '-', label="Divisive Normalization Model")
plt.legend()
plt.show()
In [23]:
list_control_observed = []
list_gabazine_observed = []
list_control_expected = []
list_gabazine_expected = []
if len(gabazine_observed_average) and len(control_observed_average):
for key in list(set(gabazine_observed_average).intersection(set(control_observed_average))):
element1, element2 = gabazine_observed_average[key], control_observed_average[key]
if not (element1 <0 or np.isclose(element1, 0, atol=tolerance) or element2<0 or np.isclose(element2, 0, atol=tolerance)):
list_gabazine_observed.append(element1)
list_control_observed.append(element2)
In [25]:
print ("GABAzine and control")
#X = np.array(list_gabazine_expected)
#y = np.array(list_gabazine_observed)
y = np.array(list_control_observed)
X = np.array(list_gabazine_observed)
idx = np.argsort(X)
X = X[idx]
y = y[idx]
linear_Model = lmfit.Model(linearModel)
DN_Model = lmfit.Model(DN_model)
lin_pars = linear_Model.make_params()
lin_result = linear_Model.fit(y, lin_pars, x=X)
lin_aic.append(lin_result.aic)
lin_chi.append(lin_result.redchi)
DN_pars = DN_Model.make_params()
DN_result = DN_Model.fit(y, DN_pars, x=X)
dn_aic.append(DN_result.aic)
dn_chi.append(DN_result.redchi)
print (lin_result.fit_report())
print (DN_result.fit_report())
ax = plt.subplot(111)
ax.scatter(X, y, alpha=0.2)
ax.set_xlim(xmin=0.)
ax.set_ylim(ymin=0.)
ax.set_xlabel("Expected")
ax.set_ylabel("Observed")
ax.set_title("Divisive Normalization and Inhibition fits")
ax.plot(X, lin_result.best_fit, '-', label="Divisive Inhibition Model")
ax.plot(X, DN_result.best_fit, '-', label="Divisive Normalization Model")
plt.show()
In [ ]: